﻿#include "FITKPlateCurvature.h"
#include "FITK_Kernel/FITKCore/FITKVec3D.h"
#include "FITK_Interface/FITKInterfaceModel/FITKElementTri.h"

namespace Interface
{
    FITKPlateCurvature::FITKPlateCurvature()
    {
        m_mesh = new FITKPlateMesh;
    }
    FITKPlateCurvature::~FITKPlateCurvature()
    {
        if (!m_mesh) delete m_mesh;
    }

    int FITKPlateCurvature::getPlateType(const QList<Core::FITKNode*>& meshPoints, const QList<Interface::FITKElementTri3*>& meshFaces)
    {
        if (!this->initPlateMesh(meshPoints, meshFaces)) return -1;

        //存储板的面积
        double plate_area = 0.0;
        //存储平板特征点面积
        double panel_area = 0.0;
        //存储单曲率板特征点面积
        double single_area = 0.0;
        //存储双曲率板特征点面积
        double double_area = 0.0;
        //计算边界点
        QSet<int> boundary = this->getBoundary();

        //遍历除边界点外的每个点，计算点的局部曲面特征
        for (OpenMesh::PolyConnectivity::VertexIter v_it = m_mesh->vertices_begin(); v_it != m_mesh->vertices_end(); ++v_it)
        {
            if (boundary.contains(v_it->idx()))continue;
            //存储各个特征点的面积
            double point_area = 0.0;
            //计算点的局部特征
            int point_feature = this->pointFeature(*v_it, point_area);
            //计算整个板的网格面积
            plate_area += point_area;
            //计算平板特征点面积
            if (point_feature == 0) panel_area += point_area;
            //计算单曲率板特征点面积
            else if (point_feature == 1) single_area += point_area;
            //计算单曲率板特征点面积
            else if (point_feature == 2) double_area += point_area;
        }

        //平面特征占比
        double panel_percentage = panel_area / plate_area;
        //单曲率面特征占比
        double single_percentage = single_area / plate_area;
        //双曲率面特征占比
        double double_percentage = double_area / plate_area;  

        //平板
        if (panel_percentage >= 0.7) 
            return 0;
        //单曲率板
        else if (single_percentage > double_percentage) 
            return 1;
        //双曲率板
        else if (double_percentage > single_percentage)
            return 2;
        //无法判定
        else return 0;
    }

    QList<double> FITKPlateCurvature::getCurvatureRadius(const int& plateType)
    {
        if (m_principalCurvature.size() == 0 || plateType == -1) return QList<double>();

        QList<double> curvature_radius{};
        //记录两个主曲率
        double k1 = m_principalCurvature[0].first;
        double k2 = m_principalCurvature[0].second;
        for (int i = 0; i < m_principalCurvature.size(); i++)
        {
            if (plateType == 1)
            {
                if ((k1 == 0 && k2 != 0) && k2 < m_principalCurvature[i].second) k2 = m_principalCurvature[i].second;
                if ((k1 != 0 && k2 == 0) && k1 < m_principalCurvature[i].first) k1 = m_principalCurvature[i].first;
            }
            if (plateType == 2)
            {
                if (k2 < m_principalCurvature[i].second) k2 = m_principalCurvature[i].second;
                if (k1 < m_principalCurvature[i].first) k1 = m_principalCurvature[i].first;
            }
        }

        if (k1 != 0 && k2 != 0)
        {
            curvature_radius.append(1 / k1);
            curvature_radius.append(1 / k2);
        }
        if (k1 != 0 && k2 == 0) curvature_radius.append(1 / k1);
        if (k1 == 0 && k2 != 0) curvature_radius.append(1 / k2);

        return curvature_radius;
    }

    bool FITKPlateCurvature::initPlateMesh(const QList<Core::FITKNode*>& meshPoints, const QList<Interface::FITKElementTri3*>& meshFaces)
    {
        if (meshPoints.empty() || meshFaces.empty()) return false;

        //对网格对象的数据清理
        m_mesh->clear();
        //记录网格点的容器
        std::vector<FITKPlateMesh::VertexHandle> vhandles{};
        //记录面顶点的容器
        std::vector<FITKPlateMesh::VertexHandle> face_vhandles{};
        //记录网格点的id
        QList<int> node_id{};

        //添加点
        for (QList<Core::FITKNode*>::const_iterator it = meshPoints.begin(); it != meshPoints.end(); it++)
        {
            FITKPlateMesh::Point m_point((*it)->x(), (*it)->y(), (*it)->z());
            //存储网格点的id
            node_id.append((*it)->getNodeID());
            vhandles.push_back(m_mesh->add_vertex(m_point));
        }

        //添加面
        for (QList<Interface::FITKElementTri3*>::const_iterator it = meshFaces.begin(); it != meshFaces.end(); ++it)
        {
            face_vhandles.clear();
            //判断是否是三角形拓扑单元类型
            if ((*it)->getEleType() != FITKModelEnum::FITKEleType::Tri3) return false;

            //获取点的id,加入到面顶点的容器中
            for (int i = 0; i < (*it)->getNodeCount(); i++)
            {
                //根据点的ID获取所在的索引
                int node_index = node_id.indexOf((*it)->getNodeID(i));
                if (node_index < 0 || node_index >= vhandles.size()) continue;
                face_vhandles.push_back(vhandles[node_index]);
            }
            //添加面
            if (face_vhandles.size() != 3) continue;

            m_mesh->add_face(face_vhandles);
        }
        return true;
    }

    QSet<int> FITKPlateCurvature::getBoundary()
    {
        QSet<int> boundary{};
        for (FITKPlateMesh::EdgeIter edge_iter = m_mesh->edges_begin(); edge_iter != m_mesh->edges_end(); ++edge_iter)
        {
            // 如果某条边是边界边
            if (m_mesh->is_boundary(*edge_iter))
            {
                //获取边的两个顶点
                FITKPlateMesh::VertexHandle vh1 = m_mesh->to_vertex_handle(m_mesh->halfedge_handle(*edge_iter, 0));
                FITKPlateMesh::VertexHandle vh2 = m_mesh->to_vertex_handle(m_mesh->halfedge_handle(*edge_iter, 1));

                //存储边的顶点索引
                boundary.insert(vh1.idx());
                boundary.insert(vh2.idx());
            }
        }
        return boundary;
    }

    int FITKPlateCurvature::pointFeature(const FITKPlateMesh::VertexHandle& vertex, double& pArea)
    {
        //记录点与相邻点角度余弦值
        QList<double> costheta{};
        //记录平均曲率算子的一部分
        OpenMesh::Vec3d mean_curvature_operator(0.0, 0.0, 0.0);
        //计算Voronoi区域面积
        double v_area = this->calVoronoiArea(vertex, costheta, mean_curvature_operator, pArea);
        //计算高斯曲率
        double gaussian_curvature = this->calGaussianCurvature(v_area, costheta);
        //计算平均曲率
        double mean_curvature = this->calMeanCurvature(v_area, mean_curvature_operator);

        //计算主曲率
        QPair<double, double> principal_curvature = this->calPrincipalCurvature(gaussian_curvature, mean_curvature);
        
        m_principalCurvature.append(principal_curvature);
        //K!=0
        if (gaussian_curvature != 0) return 2;
        //K=0
        else
        {
            //K = 0 && H == 0平面
            if (mean_curvature == 0) return 0;
            //K = 0 && H != 0抛物面
            else return 1;
        }
    }

    double FITKPlateCurvature::calTriArea(const double& a, const double& b, const double& c)
    {
        //p = 1/2(a+b+c)
        double p = 0.5*(a + b + c);
        //S = √(p(p - a)(p - b)(p - c))
        return sqrt(p*(p - a)*(p - b)*(p - c));
    }

    double FITKPlateCurvature::calVoronoiArea(const FITKPlateMesh::VertexHandle& vertex, QList<double>& costheta, OpenMesh::Vec3d& meancuroperator, double& pArea)
    {
        //Voronoi区域面积
        double v_area = 0.0;
        //获取中心点坐标
        OpenMesh::Vec3d center_point = m_mesh->point(vertex);

        //获取相邻点,遍历包含中心点的所有面
        for (FITKPlateMesh::VertexFaceIter vf_it = m_mesh->vf_iter(vertex); vf_it.is_valid(); ++vf_it)
        {
            //记录相邻点
            int i = 0;
            OpenMesh::Vec3d neighbor_point[2];
            neighbor_point[0] = OpenMesh::Vec3d(0, 0, 0);
            neighbor_point[1] = OpenMesh::Vec3d(0, 0, 0);

            //遍历面的所有顶点
            for (FITKPlateMesh::FaceVertexIter fv_it = m_mesh->fv_iter(*vf_it); fv_it.is_valid(); ++fv_it)
            {
                if (*fv_it != vertex)
                {
                    neighbor_point[i] = m_mesh->point(*fv_it);
                    i++;
                }
            }
            //计算三角形三个角度余弦值OPQ
            double lengthp = (neighbor_point[0] - center_point).length();
            double lengthq = (neighbor_point[1] - center_point).length();
            double lengtho = (neighbor_point[1] - neighbor_point[0]).length();
            //cos o = P·Q/|P||Q|
            double dotpq = (neighbor_point[0] - center_point) | (neighbor_point[1] - center_point);
            double coso = dotpq / (lengthp * lengthq); 

            //cos p = O·Q/|O||Q|
            double dotoq = (center_point - neighbor_point[1]) | (neighbor_point[0] - neighbor_point[1]);
            double cosp = dotoq / (lengtho * lengthq);
            //cos q = O·P/|O||P|
            double dotop = (neighbor_point[1] - neighbor_point[0]) | (center_point - neighbor_point[0]);
            double cosq = dotop / (lengtho * lengthp);

            //if (coso == 0 || cosp == 0 || cosq == 0) continue;

            //将角度余弦值添加到容器中
            costheta.append(coso);

            //计算p、q的余切值
            double cotp = cosp / sqrt(1 - cosp * cosp);
            double cotq = cosq / sqrt(1 - cosq * cosq);

            //平均曲率计算算子（OQ*cot∠P + OP*cot∠Q）
            meancuroperator += (neighbor_point[0] - center_point) * cotp + (neighbor_point[1] - center_point) * cotq;

            //三角形面积
            pArea += calTriArea(lengtho, lengthp, lengthq);

            //如果存在钝角或相邻点角为直角，则Voronoi不规则
            if (coso < 0 || cosp <= 0 || cosq <= 0)
            {
                //如果中心点角度为钝角，Voronoi区域面积为1/2三角形面积
                if (coso < 0)
                {
                    double area_T = 0.5 * calTriArea(lengtho, lengthp, lengthq);
                    v_area += area_T;
                }
                //如果其他点角度为钝角，Voronoi区域面积为1/4三角形面积
                else
                {
                    double area_T = 0.25 * calTriArea(lengtho, lengthp, lengthq);
                    v_area += area_T;
                }
            }
            else
            {
                //一个三角形内Voronoi区域面积S =1/8（|OQ|²cot∠P + |OP|²cot∠Q）
                double area_T = 0.125 * (lengthp * lengthp * cotp + lengthq * lengthq * cotq);
                v_area += area_T;
            }
        }
        //Voronoi面积
        return v_area;
    }

    double FITKPlateCurvature::calGaussianCurvature(const double& vArea, const QList<double>& costheta)
    {
        //计算中心点角度和
        double o_theta = 0.0;
        for (int i = 0; i < costheta.size(); i++)
        {
            o_theta += acos(costheta[i]);
        }

        //计算2pi-∑theta,设置误差,如果角度小于0.01°则认定为0°
        double angela = 2 * M_PI - o_theta;
        if (angela < 1.7453292519943296e-4) angela = 0;

        return(angela / vArea);
    }

    double FITKPlateCurvature::calMeanCurvature(const double& vArea,OpenMesh::Vec3d& meanCurOperator)
    {
        if (abs(meanCurOperator[0]) < 1e-3) meanCurOperator[0] = 0.0;
        if (abs(meanCurOperator[1]) < 1e-3) meanCurOperator[1] = 0.0;
        if (abs(meanCurOperator[2]) < 1e-3) meanCurOperator[2] = 0.0;

        //计算平均曲率算子 K(xi) = 0.5 * meanCurOperator / vArea
        OpenMesh::Vec3d Ki = 0.5 * (meanCurOperator / vArea);
        //计算平均曲率 H = 0.5 * ||K(xi)||
        return 0.5 * Ki.norm();
    }

    QPair<double, double> FITKPlateCurvature::calPrincipalCurvature(const double& gaussianCurvature, const double& meanCurvature)
    {
        double k_operator = meanCurvature * meanCurvature - gaussianCurvature;

        if (k_operator < 0) return QPair<double, double>(0, 0);

        //k1 = H + √(H²-K)
        double k1 = meanCurvature + sqrt(k_operator);
        //k2 = H - √(H²-K)
        double k2 = meanCurvature - sqrt(k_operator);

        return QPair<double, double>(k1, k2);
    }
}